
clear Mcomp3  CL CR Lnexp_cumulative_IT     LnexpXX  Lnexp_cumulative MeltType Ln_cumulative


if Gar2SpTransition~=Pb(Gar2SpTransition_index)
    'Gar2Sp transition pressure mismatch'
    Gar2SpTransition_index=size(Ts,2)-1;
else
    
    
end

if Sp2PlagTransition~=Pb(Sp2PlagTransition_index)
    Sp2PlagTransition_index=size(Ts,2)-1;
    'Sp2Plag transition pressure mismatch'
end


if Gar2SpTransition_index==size(Ts,2)
    Tdiff = NaN;
else
    Tdiff =  Ts(Gar2SpTransition_index) - Ts(Gar2SpTransition_index+1);
end

if porosity == 1
    Lnexp = ln;
    BAI4 = BAI3(:,1:9);
    CL_OUT = CL_BAI3;
end



CR_OUT=[Petrogen_Co';  CR_OUT(1:end-1,:)];

BulkP = P_out(:,U_index)./P_out(:,Th_index);
BulkDo = Do_out(:,U_index)./Do_out(:,Th_index);
BulkD = Dn_Test(:,U_index)./Dn_Test(:,Th_index);
BulkD_Ra = Dn_Test(:,Th_index)./Dn_Test(:,Ra_index);
Ds2save = [DUTh_out BulkD];


%%
ThU_Source = Petrogen_Co(Th_index)./Petrogen_Co(U_index);
RaTh_Source = Petrogen_Co(Ra_index)./Petrogen_Co(Th_index);

ZeroIngrowth_NF_Th230 = (CL_OUT(:,Th_index)./CL_OUT(:,U_index))./(ThU_Source);
ZeroIngrowth_NF_Ra226 = (CL_OUT(:,Ra_index)./CL_OUT(:,Th_index))./(RaTh_Source);

CompleteIngrowth_NF_Th230  = (CL_OUT(:,Th_index)./CL_OUT(:,U_index))./(CR_OUT(:,Th_index)./CR_OUT(:,U_index));
CompleteIngrowth_NF_Ra226= (CL_OUT(:,Ra_index)./CL_OUT(:,Th_index))./(CR_OUT(:,Ra_index)./CR_OUT(:,Th_index));

%%
% upwelling rate = 2/pi * full spreading rate
% t=25.9e5;	%0.5cm/yr
% t=25.9e3;	%0.5cm/yr
% t=13.0e3;	%1cm/yr
%  t=6.5e3;	%2cm/yr
% t=4.3e3;	%3cm/yr
%  t=3.2e3;	%4cm/yr
%  t=2.6e3;	%5cm/yr
% t=2.2e3;	%6cm/yr
%  t=1.9e3;	%7cm/yr
% t=1.6e3;	%8cm/yr
%  t=1.4e3;	%9cm/yr
% t=1.3e3;	%10cm/yr
% halflifeTh230 = 7.56e4;
% lambdaTh230 = log(2)./halflifeTh230;
%
%
% Residues = (CR_OUT(:,Th_index)./CR_OUT(:,U_index));
% Melts = (CL_OUT(:,Th_index)./CL_OUT(:,U_index));
% Residues_shift = Residues;
% Residues_shift=[NaN; Residues(1:end-1)];
%
%
% testtime = 1:(size(CompleteIngrowth_NF_Th230,1));
% TimeSinceMelting=CompleteIngrowth_NF_Th230;
% TimeSinceMelting(~isnan(CompleteIngrowth_NF_Th230))=1;
% testtime = [zeros(1,TimeSinceMelting-1) 1:(size(Pb,2)-TimeSinceMelting+1)];
% testtime=t.*testtime;
% IntermediateUpwellingSource=ThU_Source+((Residues-ThU_Source).*(1-exp(-1*lambdaTh230.*testtime')));
%
%
%
% TimeSinceMelting = find(~isnan(CompleteIngrowth_NF_Th230), 1);
%
% IngrownMantle(1)=ThU_Source;
%
% for yPb = 2:size(Pb)
%     IngrownMantle(yPb)=3;
% end
%
%
%
% Fraction= (1-exp(-1*lambdaTh230.*testtime'));
% Fraction = 0.5;
% IntermediateUpwellingSource = ThU_Source.*(1-Fraction)+(Residues.*Fraction);
% IntermediateUpwellingSource_Excess =  (CL_OUT(:,Th_index)./CL_OUT(:,U_index))./IntermediateUpwellingSource;
%
%
% % figure()
% % hold on
% %
% % plot(1./Melts, 1./ThU_Source.*ones(size(Melts,1),1),'k')
% % plot(1./Residues, 1./ThU_Source.*ones(size(Melts,1),1),'r')
% % plot(1./Melts(TimeSinceMelting), 1./ThU_Source,'ko')
% % plot(1./Residues(TimeSinceMelting), 1./ThU_Source,'ro')
% %
% % plot(1./Melts, 1./Residues.*ones(size(Melts,1),1),'g')
% % %plot(1./Residues, 1./ThU_Source.*ones(size(Melts,1),1),'r')
% % plot([.2 .75],[.2 .75],'k-')
% % xlabel('U/Th')
% % ylabel('230Th/232Th')
% % axis square
% % box on
%
% figure()
% hold on
% plot(CR_OUT(:,Th_index)./CR_OUT(:,U_index),Pb, 'ko-')
% %plot(CL_OUT(:,Th_index)./CL_OUT(:,U_index),Pb, 'go-')
%
% % plot(ZeroIngrowth_NF_Th230,Pb, 'bo-')
% % plot(CompleteIngrowth_NF_Th230,Pb, 'Petrogen_Co-')
%
% plot(IntermediateUpwellingSource,Pb, 'yo-')
%
% vline(ThU_Source,'r-')
% vline(1,'k-')
%
%
% set(gca,'YDir','Reverse')
%
%
%
% figure()
% hold on
% plot(CR_OUT(:,Th_index),Pb, 'ko-')
% plot(CR_OUT(:,U_index),Pb, 'ro-')
%
%
% xline(Petrogen_Co(Th_index),'k-')
% xline(Petrogen_Co(U_index),'r-')
% legend('Th','U','Th S','U S')
% set(gca,'YDir','Reverse')
%
%
% figure()
% hold on
% plot(ZeroIngrowth_NF_Th230,Pb, 'bo-')
% plot(CompleteIngrowth_NF_Th230,Pb, 'Petrogen_Co-')
% plot(IntermediateUpwellingSource_Excess,Pb, 'yo-')
%
% vline(1,'k-')
%
% legend('zero','complete','intermediate')
%
% set(gca,'YDir','Reverse')
%
%
% return

%%
for z = 1:size(Lnexp,2)
    meltsSoFar = Lnexp(1:z)./sum(Lnexp(1:z));
    AvgTh_z(z) = CL_OUT(1:z,Th_index)'*meltsSoFar';
    AvgU_z(z) = CL_OUT(1:z,U_index)'*meltsSoFar';
    AvgRa_z(z) = CL_OUT(1:z,Ra_index)'*meltsSoFar';
    AvgTh_R_z(z) = CR_OUT(1:z,Th_index)'*meltsSoFar';
    AvgU_R_z(z) = CR_OUT(1:z,U_index)'*meltsSoFar';
    AvgRa_R_z(z) = CR_OUT(1:z,Ra_index)'*meltsSoFar';
end

ZeroIngrowth_Pool_Th230 = (AvgTh_z./AvgU_z)./ThU_Source;
ZeroIngrowth_Pool_Ra226 = (AvgRa_z./AvgTh_z)./RaTh_Source;
CompleteIngrowth_Pool_Th230 = (AvgTh_z./AvgU_z)./(AvgTh_R_z./AvgU_R_z);
CompleteIngrowth_Pool_Ra226= (AvgRa_z./AvgTh_z)./(AvgRa_R_z./AvgTh_R_z);

DISEQ_meltmass = [ZeroIngrowth_Pool_Th230(end) ZeroIngrowth_Pool_Ra226(end) CompleteIngrowth_Pool_Th230(end) CompleteIngrowth_Pool_Ra226(end)];

%% Figures to check out diseq
% figure()
% hold on
% box on
% goodmelts = find(Lnexp>0);
% yline(1,'HandleVisibility','off');
%
% plot(CL_OUT(goodmelts,Th_index), ZeroIngrowth_NF_Th230(goodmelts),'r-','LineWidth',2)
% plot(CL_OUT(goodmelts,Th_index), ZeroIngrowth_Pool_Th230(goodmelts),'r:','LineWidth',3)
%
% plot(CL_OUT(goodmelts,Th_index), CompleteIngrowth_NF_Th230(goodmelts),'k-','LineWidth',2)
% plot(CL_OUT(goodmelts,Th_index), CompleteIngrowth_Pool_Th230(goodmelts),'k:','LineWidth',3)
%
% legend('Zero Ingrowth NF', 'Zero Ingrowth Pool' ,'Complete Ingrowth NF', 'Complete Ingrowth Pool','Location','Best')
% xlabel('Th ppm')
% ylabel('(230Th)/(238U)')
%
%
% figure()
% hold on
% box on
% goodmelts = find(Lnexp>0);
% yline(1,'HandleVisibility','off');
%
% plot(CL_OUT(goodmelts,Th_index), ZeroIngrowth_NF_Ra226(goodmelts),'r-','LineWidth',2)
% plot(CL_OUT(goodmelts,Th_index), ZeroIngrowth_Pool_Ra226(goodmelts),'r:','LineWidth',3)
%
% plot(CL_OUT(goodmelts,Th_index), CompleteIngrowth_NF_Ra226(goodmelts),'k-','LineWidth',2)
% plot(CL_OUT(goodmelts,Th_index), CompleteIngrowth_Pool_Ra226(goodmelts),'k:','LineWidth',3)
%
% legend('Zero Ingrowth NF', 'Zero Ingrowth Pool' ,'Complete Ingrowth NF', 'Complete Ingrowth Pool','Location','Best')
% xlabel('Th ppm')
% ylabel('(226Ra)/(230Th)')

%% Figures to check out Ds
% figure()
% hold on
% box on
% plot(CumFrac, Do_out(:,Th_index),'k.-')
% plot(CumFrac, Do_out(:,U_index),'r.-')
% plot(CumFrac,DUTh_out(:,2),'k:');
% plot(CumFrac,DUTh_out(:,1),'r:');
% legend('Bulk DTh','Bulk DU','CPX DTh','CPX DU')
% xlabel('Near-Frac F')
% ylabel('D')
%
% figure()
% hold on
% box on
% plot(CumFrac, CompleteIngrowthIncrementalMelts,'k.-')
% plot(CumFrac, BulkD,'r.-')
% legend('Complete Ingrowth 230Th','Bulk DU/DTh')
% xlabel('Near-Frac F')
%
% figure()
% hold on
% box on
% plot(CumFrac, CompleteIngrowthIncrementalMelts_226Ra,'k.-')
% plot(CumFrac, BulkD_Ra,'r.-')
% legend('Complete Ingrowth 226Ra','Bulk DTh/DRa')
% xlabel('Near-Frac F')


%%

% Outputs:
% XM_OUT - Mantle mode [Aug, Opx, Oliv, Plag, Sp]
% XOW_OUT - Mantle composition [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% BAI3 - All near-fractional melts [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% BAI4 - Extracted near-fractional melts [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% Ts - Solidus temperature (C)
% Tsdry - Dry solidus temperature (C)
% T - Temperature corrected for latent heat of melting based on the
%       amount of melt extracted at each step (C)
% CumFrac, IncFrac - cumulative and near-fractional melt fraction
% Fc  - max(CumFrac); includes melt contained in residue
% rXOW - Residual mantle composition (residue = retained melt + solid) [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% MgRatio -  [Mg# residue, Mg# melt] ;
% NaKratio - NaK# melt
% test - [Pb Ts Tb (1-XMGliq) WNA WT WK WNACA Sp2PlagTransitionPressure Gar2SpTransitionPressure QCA];
% initialBulkComp - the mantle bulk composition used
% CL_OUT - trace element composition of extracted near-fractional melts
% CS_OUT - trace element composition of the solid
% CR_OUT - trace element composition of the residue
%           (residue = retained melt + solid)
% CL_BAI3 - trace element composition of all near-fractional melts
% CS_minerals_OUT - trace element composition minerals in solid residue, right now it is just CPX
% Do_out - the bulk modal-melting D at each step
% P_out - the non-modal melting P at each step
% Dn_Test - the non-modal melting D at each step
% ln - mass of new liquid generated at each step
% WnS - mass of the solid in the residue at each step
% Ln - mass of the liquid at each step
% Qn - fraction of the melt that remains in the residue at each step
% Lnr - mass of melt that remains in the residue  at each step
% Lnexp - mass of melt expelled at each step
% Wnr - mass of residue at each step
% Gar2SpTransition_index - index of the garnet to spinel phase transition
%                           (i.e, phase transition pressure = Pb(Gar2SpTransition_index)
% Sp2PlagTransition_index - index of the spinel to plagioclase phase transition
%                           (i.e, phase transition pressure = Pb(Sp2PlagTransition_index)
% Lnexp_IT -  mass of melt expelled in 2D
% DUTh_out -  the bulk DU/DTh at each step (which is pressure dependent)

assignin('caller', sprintf('XM_OUT_%g', runnumber),XM_OUT)
assignin('caller', sprintf('XOW_OUT_%g', runnumber), XOW_OUT)
assignin('caller', sprintf('BAI3_%g', runnumber),BAI3)
assignin('caller', sprintf('BAI4_%g', runnumber),BAI4)
assignin('caller', sprintf('Ts_%g', runnumber),Ts)
assignin('caller', sprintf('T_%g', runnumber), T)
assignin('caller', sprintf('Tb_%g', runnumber),Tb)
assignin('caller', sprintf('CumFrac_%g', runnumber),CumFrac)
assignin('caller', sprintf('IncFrac_%g', runnumber), IncFrac)
assignin('caller', sprintf('Fc_%g', runnumber),Fc)
assignin('caller', sprintf('rXOW_%g', runnumber),rXOW)
assignin('caller', sprintf('MgRatio_%g', runnumber), MgRatio)
assignin('caller', sprintf('NaKratio_%g', runnumber),NaKratio)
assignin('caller', sprintf('test_%g', runnumber),test)
assignin('caller', sprintf('initialBulkComp_%g', runnumber), initialBulkComp)
assignin('caller', sprintf('CL_OUT_%g', runnumber),CL_OUT)
assignin('caller', sprintf('CS_OUT_%g', runnumber), CS_OUT)
assignin('caller', sprintf('CR_OUT_%g', runnumber), CR_OUT)
assignin('caller', sprintf('CL_BAI3_%g', runnumber),CL_BAI3)
assignin('caller', sprintf('CS_minerals_OUT_%g', runnumber),CS_minerals_OUT)
assignin('caller', sprintf('Do_out_%g', runnumber),Do_out)
assignin('caller', sprintf('P_out_%g', runnumber),P_out)
assignin('caller', sprintf('Dn_Test_%g', runnumber), Dn_Test)
assignin('caller', sprintf('ln_%g', runnumber), ln)
assignin('caller', sprintf('WnS_%g', runnumber), WnS)
assignin('caller', sprintf('Ln_%g', runnumber), Ln)
assignin('caller', sprintf('Qn_%g', runnumber), Qn)
assignin('caller', sprintf('Lnr_%g', runnumber), Lnr)
assignin('caller', sprintf('Lnexp_%g', runnumber), Lnexp)
assignin('caller', sprintf('Wnr_%g', runnumber),Wnr)
assignin('caller', sprintf('Gar2SpTransition_index_%g', runnumber), Gar2SpTransition_index)
assignin('caller', sprintf('Sp2PlagTransition_index_%g', runnumber), Sp2PlagTransition_index)






%%

%for 2D melting calculation
%assumes the phase transitions occur at the same pressure all across, which is a
%simiplifcation. In reality, at the wings, off-axis, garnet will be stable to more shallow
%pressures
secyr = 3600*24*365.25;
M = Lnexp_IT;
Mr = Iv_2D'.*(M'/dz);  % Melt production rate c(1/s)
Mr = Mr';
% Calculate melt productivity and crustal thickness
Mr(Mr<0) = 0;      % Eliminate values where melt productivity is negative
tMx = sum(Mr).*dz;  % total melt at each X-value
% crustal thickness changed 7/7/17 to divide by full spreading rate, not half spreading rate
ct = sum(tMx).*dx./(10*(2.*afewfunrnuns_U_vary(afewfunrnuns_U)/secyr));


[minValue,closestIndex] = min(abs(T-600));
BrittleDuctileTransPress = Pb(closestIndex);

AAdydx = diff(T)./diff(Pb);
[cq,dq]=find(AAdydx>6);
[aq,bq]=(max(T(dq)-Ts(dq),[],2));
LastAdiabaticMeltPressure = Pb(dq(bq));
if isempty(LastAdiabaticMeltPressure)==1
    LastAdiabaticMeltPressure=nan;
end

[minValue,closestIndex] = min(abs(AAdydx-20));
LastAdiabaticMeltPressure_v2 = Pb(closestIndex);


[minValue,closestIndex] = min(abs(T-1200));
LastAdiabaticMeltPressure_v3 = Pb(closestIndex);

AdiabaticMeltingInfo = [LastAdiabaticMeltPressure LastAdiabaticMeltPressure_v2 LastAdiabaticMeltPressure_v3 BrittleDuctileTransPress];

%%
%         close all
%
%         figure(); hold on;
%         plot(T,Pb,'b');
%         plot(Ts,Pb,'r');
%         yline(LastAdiabaticMeltPressure_v2,'k-','dTdP=20C/kbar');
%         yline(LastAdiabaticMeltPressure,'k-','max T-Ts')
%         yline(LastAdiabaticMeltPressure_v3,'k-','1200C isotherm')
%         set(gca,'YDir','Reverse')
%
%
%         figure(); hold on;
%         plot(3.154e9.*Iv_2D(:,101),Pb,'k-');
%         yline(LastAdiabaticMeltPressure);
%         yline(LastAdiabaticMeltPressure_v2);set(gca,'YDir','Reverse')
%
%         yline(LastAdiabaticMeltPressure_v2,'k-','dTdP=20C/kbar');
%         yline(LastAdiabaticMeltPressure,'k-','max T-Ts')
%         yline(LastAdiabaticMeltPressure_v3,'k-','1200C isotherm')
%         set(gca,'YDir','Reverse')
%
%
%
%
%
%         tempdiff = Iv_2D.*3.154e9;
%         %tempdiff(tempdiff<0) = NaN;
%
%         upwellingThres = 1;%0.04;
%         [minValue,closestIndex] = min(abs(tempdiff-upwellingThres));
%         PressureSlow = Pb(closestIndex);
%
%
%         [minValue,closestIndex] = min(abs(IT_extended-1200));
%         LAB_1200_Iso= Pb(closestIndex);
%
%
%
%         for kklm= 1:size(IT_extended,2)
%             T = IT_extended(:,kklm)';
%             AAdydx = diff(T)./diff(Pb);
%             [cq,dq]=find(AAdydx>6);
%             [aq,bq]=(max(T(dq)-Ts(dq),[],2));
%             LastAdiabaticMeltPressure_2D(kklm) = Pb(dq(bq));
%
% %            LastMeltPressure_2D(kklm) = Pb(max(find(Lnexp_IT(:,kklm)>0)));%Pb(find(Lnexp_IT(:,kklm)==0,1,'last'));
%
%         end
%
%
%         figure()
%         pcolor(x/1e3,Pb',tempdiff);
%         shading flat
%        % caxis manual
%         %caxis([0 160]);
%         colorbar
%         colormap('jet')
%         hold on
%         %plot(x/1e3,Gar2SpTransition.*ones(size(x)),'r--', 'LineWidth',3)
%         plot(x/1e3,Sp2PlagTransition.*ones(size(x)),'k-', 'LineWidth',3)
%         %plot(x/1e3,Pbextractedindex,'Color','r', 'LineWidth',3)
%
%         plot(x/1e3,PressureSlow,'c-', 'LineWidth',3)
%         plot(x/1e3,LAB_1200_Iso,'y-', 'LineWidth',3)
%
%         plot(x/1e3,LastAdiabaticMeltPressure_2D,'r:', 'LineWidth',3)
%
% %        plot(x/1e3,LastMeltPressure_2D,'r-', 'LineWidth',3)
%
%
%         legend('velocity','Sp2Plag','Local-Element-Frac','1200C Isotherm','Last Adiabatic','LastMelt')
%
%         vline(x(101-50)/1e3)
%         vline(x(101-90)/1e3)
%
%         vline(x(101+50)/1e3)
%         vline(x(101+90)/1e3)
%         ylabel('Pressure (kbars)')
%         xlabel('Distance from axis (km)')
%            axis([-200 200 0 15])
%         box on
%         set(gca,'XMinorTick','on','YMinorTick','on')
%         set(gca, 'ydir', 'reverse');
%
%
%
%
%
%         figure()
%         tempdiff=Mr;
%         tempdiff(tempdiff==0) = NaN;
%         sdp=pcolor(x/1e3,Pb',tempdiff);
%         alpha(sdp,0.6);
%         shading flat
%        % caxis manual
%         %caxis([0 160]);
%         colorbar
%         colormap('jet')
%         hold on
%         %plot(x/1e3,Gar2SpTransition.*ones(size(x)),'r--', 'LineWidth',3)
%         plot(x/1e3,Sp2PlagTransition.*ones(size(x)),'k-', 'LineWidth',3)
%         plot(x/1e3,ct/3.08.*ones(size(x)),'k-', 'LineWidth',3,'Color',grey5)
%         %plot(x/1e3,Pbextractedindex,'Color','r', 'LineWidth',3)
%
%         plot(x/1e3,PressureSlow,'c-', 'LineWidth',3,'Color',black)
%         plot(x/1e3,LAB_1200_Iso,'y-', 'LineWidth',3)
%
%         plot(x/1e3,LastAdiabaticMeltPressure_2D,'r-', 'LineWidth',3)
%         %plot(x/1e3,LastMeltPressure_2D,'r-', 'LineWidth',3)
%
%
%         legend('melt production rate','Sp2Plag','crust','Local-Element-Frac','1200C Isotherm','Last Adiabatic')
%
%         vline(x(101-50)/1e3)
%         vline(x(101-90)/1e3)
%
%         vline(x(101+50)/1e3)
%         vline(x(101+90)/1e3)
%         ylabel('Pressure (kbars)')
%         xlabel('Distance from axis (km)')
%            axis([-200 200 0 15])
%         box on
%         set(gca,'XMinorTick','on','YMinorTick','on')
%         set(gca, 'ydir', 'reverse');
%
%
%
%         figure(); hold on; plot(Iv_2D(:,101),Pb,'k-'); yline(LastAdiabaticMeltPressure); yline(LastAdiabaticMeltPressure_v2);set(gca,'YDir','Reverse')
%
%   return
%         figure(); hold on; plot(AAdydx,Pb(2:end),'k-'); yline(LastAdiabaticMeltPressure); yline(LastAdiabaticMeltPressure_v2);set(gca,'YDir','Reverse')
%



%%


% Calculate composition of pooled melts
% Total pooled melt from all incremental melts produced in model space

for aa = 1:201
    Mcomp3(aa,:,:) = BAI4;
    CL(aa,:,:) = CL_OUT;
    CR(aa,:,:) = CR_OUT;
end

for i = 1:9
    pool_full(i) = sum(sum(Mr.*Mcomp3(:,:,i)'))/sum(Mr(:));
end

for i = 1:size(PetrogenTraceElement_Strings,2)
    poolCL_full(i) = sum(sum(Mr.*CL(:,:,i)'))/sum(Mr(:));
    poolCR_full(i) = sum(sum(Mr.*CR(:,:,i)'))/sum(Mr(:));
end


ZeroIngrowth_Pool_Th230_full =       (poolCL_full(Th_index)./poolCL_full(U_index))./ThU_Source;
ZeroIngrowth_Pool_Ra226_full  =      (poolCL_full(Ra_index)./poolCL_full(Th_index))./RaTh_Source;
CompleteIngrowth_Pool_Th230_full  =  (poolCL_full(Th_index)./poolCL_full(U_index))./(poolCR_full(Th_index)./poolCR_full(U_index));
CompleteIngrowth_Pool_Ra226_full  =  (poolCL_full(Ra_index)./poolCL_full(Th_index))./(poolCR_full(Ra_index)./poolCR_full(Th_index));





for lll = 1:size(Lnexp_IT,2)
    LnexpXX = Lnexp_IT(:,lll);
    for kl = 1:size(LnexpXX,1)
        Lnexp_cumulative(kl) = sum(LnexpXX(1:kl)); % using (1-Wnr) has floating point issues, so use this instead
    end
    Lnexp_cumulative(LnexpXX==0)=0;
    Lnexp_cumulative_IT(:,lll) = Lnexp_cumulative;
end
FP_Full =  sum(sum(Mr.*Lnexp_cumulative_IT))/sum(Mr(:));



%COLUMN POOLED MELTS
for ll=1:size(Lnexp_cumulative_IT,2)
    variablypooledLnExp(ll) =  (nansum(Mr(:,ll).*Lnexp_cumulative_IT(:,ll)))/nansum(Mr(:,ll));
end

for ll=1:size(Lnexp_cumulative_IT,2)
    variablypooledPb(ll) =  (nansum(Mr(:,ll).*Pb'))/nansum(Mr(:,ll)); %summing each column
end

for ll=1:size(Lnexp_cumulative_IT,2)
    variablypooledTs(ll) =  (nansum(Mr(:,ll).*Ts'))/nansum(Mr(:,ll)); %summing each column
    
    columnME = squeeze(Mcomp3(ll,:,:));
    Pool_column_ME(ll,:)=(nansum(Mr(:,ll).*columnME))/nansum(Mr(:,ll)); %summing each column
    
    columnTE = squeeze(CL(ll,:,:));
    Pool_column_TE(ll,:)=(nansum(Mr(:,ll).*columnTE))/nansum(Mr(:,ll)); %summing each column
end


columnsWithMelt = find(variablypooledLnExp>0);
columnsWithMelt_ALL=columnsWithMelt;
columnsWithMelt(columnsWithMelt>101)=[];

% Pooled melt limiting off-axis migration distance
for k = 0:1:length(x)/2
    for i = 1:9
        poolME_ininc(k+1,i) = sum(sum(Mr(:,1+k:end-k).* ...
            Mcomp3(1+k:end-k,:,i)'))/sum(sum(Mr(:,1+k:end-k)));
    end
    
    for i = 1:size(CL,3)
        poolTE_ininc(k+1,i) = sum(sum(Mr(:,1+k:end-k).* ...
            CL(1+k:end-k,:,i)'))/sum(sum(Mr(:,1+k:end-k)));
        
        poolCR_ininc(k+1,i) = sum(sum(Mr(:,1+k:end-k).* ...
            CR(1+k:end-k,:,i)'))/sum(sum(Mr(:,1+k:end-k)));
    end
    
    for i = 1
        Pb_Melting_ininc(k+1) = sum(sum(Mr(:,1+k:end-k).*Pb'))/sum(sum(Mr(:,1+k:end-k)));
        T_Melting_ininc(k+1) = sum(sum(Mr(:,1+k:end-k).*Ts'))/sum(sum(Mr(:,1+k:end-k)));
    end
    
    for i = 1
        Lnexp_cumulative_IT_ininc(k+1,i) = sum(sum(Mr(:,1+k:end-k).*Lnexp_cumulative_IT(:,1+k:end-k)))/sum(sum(Mr(:,1+k:end-k)));
    end
    
end


for ll = 1:size(poolCR_ininc,1)
    ZeroIngrowth_Pool_Th230_ininc(ll) =       (poolTE_ininc(ll,Th_index)./poolTE_ininc(ll,U_index))./ThU_Source;
    ZeroIngrowth_Pool_Ra226_ininc(ll)  =      (poolTE_ininc(ll,Ra_index)./poolTE_ininc(ll,Th_index))./RaTh_Source;
    CompleteIngrowth_Pool_Th230_ininc(ll)  =  (poolTE_ininc(ll,Th_index)./poolTE_ininc(ll,U_index))./(poolCR_ininc(ll,Th_index)./poolCR_ininc(ll,U_index));
    CompleteIngrowth_Pool_Ra226_ininc(ll)  =  (poolTE_ininc(ll,Ra_index)./poolTE_ininc(ll,Th_index))./(poolCR_ininc(ll,Ra_index)./poolCR_ininc(ll,Th_index));
end

%finds cases where garnet is not in the residue
garnetpossible = 1:Gar2SpTransition_index;
garnetexist= find(XM_OUT(garnetpossible,6)>0);

for ll = 1:size(Mr,2)
    columnSumGarnetReal_Weight(ll,:) = sum(Mr(garnetexist,ll));
    columnSumGarnet_Weight(ll,:) = sum(Mr(1:Gar2SpTransition_index,ll));
    columnSumSpinel_Weight(ll,:) = sum(Mr(Gar2SpTransition_index+1:Sp2PlagTransition_index,ll));
    columnSumPlag_Weight(ll,:) = sum(Mr(Sp2PlagTransition_index+1:end,ll));
    columnSumWeight(ll) = sum(Mr(:,ll));
end

for k = 0:1:length(x)/2
    ct2_ininc(k+1)=sum(columnSumWeight(k+1:end-k)*dz).*dx./(2.*10*afewfunrnuns_U_vary(afewfunrnuns_U)/secyr);
    
    GarnetReal_ct2_ininc(k+1)=sum(columnSumGarnetReal_Weight(k+1:end-k)*dz).*dx./(2.*10*afewfunrnuns_U_vary(afewfunrnuns_U)/secyr);
    
    Garnet_ct2_ininc(k+1)=sum(columnSumGarnet_Weight(k+1:end-k)*dz).*dx./(2.*10*afewfunrnuns_U_vary(afewfunrnuns_U)/secyr);
    Spinel_ct2_ininc(k+1)=sum(columnSumSpinel_Weight(k+1:end-k)*dz).*dx./(2.*10*afewfunrnuns_U_vary(afewfunrnuns_U)/secyr);
    Plag_ct2_ininc(k+1)=sum(columnSumPlag_Weight(k+1:end-k)*dz).*dx./(2.*10*afewfunrnuns_U_vary(afewfunrnuns_U)/secyr);
end



xkm = x./1000;
distance_ininc = abs(xkm(1:101));


allThickness = [Garnet_ct2_ininc; Spinel_ct2_ininc; Plag_ct2_ininc; ct2_ininc];
percentCrust = [Garnet_ct2_ininc./ct2_ininc; Spinel_ct2_ininc./ct2_ininc; Plag_ct2_ininc./ct2_ininc];



if XM_OUT(1,6) ==0
    'no initial garnet'
    percentCrust(:,:) = 0;
    
else if XM_OUT(Gar2SpTransition_index-1,6) ==0
        'garnet is exhausted'
        Garnet_proportion = GarnetReal_ct2_ininc./ct2_ininc;
        percentCrust(1,:) = Garnet_proportion;
        percentCrust(2:3,:) = 0;
    end
end

%%



Iv = Iv_2D(:,columnnum);

%given the PETROGEN 1D results, calculate useful post-processing values:
[ExtFrac, aluminumOutIndex, cpxOutIndex,TransitionSp2Plag, TransitionGar2Sp,aluminumOutFrac, cpxOutFrac,MeltExtractedPb, BAI3, ...
    XOW_OUT, rXOW,initialBulkComp,FirstMeltPb,FirstAt2percent,FirstMeltIndex,FirstAt2Index,LastMeltPb,LastMeltIndex,crustalthickness,...
    poolCL_Iv,poolCL,poolBAI4_Iv,poolBAI4,running_mBAI4,running_CL_OUT,runningPressure] = PETROGEN2019_post_proc(...
    Pb, eval(sprintf('T_%g', runnumber)),  eval(sprintf('ln_%g', runnumber)), eval(sprintf('Lnexp_%g', runnumber)), eval(sprintf('WnS_%g', runnumber)), eval(sprintf('Wnr_%g', runnumber)), eval(sprintf('Lnr_%g', runnumber)), ...
    eval(sprintf('XM_OUT_%g', runnumber)), eval(sprintf('BAI3_%g', runnumber)), eval(sprintf('XOW_OUT_%g', runnumber)), eval(sprintf('rXOW_%g', runnumber)),...
    eval(sprintf('initialBulkComp_%g', runnumber)),eval(sprintf('CumFrac_%g', runnumber)), eval(sprintf('Gar2SpTransition_index_%g', runnumber)), eval(sprintf('Sp2PlagTransition_index_%g', runnumber)),...
    Iv, eval(sprintf('CL_OUT_%g', runnumber)), eval(sprintf('BAI4_%g', runnumber)),...
    porosity,Petrogen_Dmatrix,eval(sprintf('CL_BAI3_%g', runnumber)));


% Output:

% ExtFrac
% aluminumOutIndex = Index at aluminous phase exhuastion
% cpxOutIndex = Index at cpx phase exhuastion
% TransitionSp2Plag= pressure of the spinel to plagioclase phase transition
% TransitionGar2Sp = pressure of the garnet to spinel phase transition
% aluminumOutFrac = Melt fraction at aluminous phase exhuastion
% cpxOutFrac= Melt fraction at cpx exhuastion
% MeltExtractedPb  = Presure of first melt extraction
% BAI3 = All incremental melts [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O,Mg#,CaO/Al2O3, NaK#]
% rXOW = Residual Mantle Bulk Composition [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O,Mg#,CaO/Al2O3, NaK#]
% initialBulkComp = Mantle Bulk Composition [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O,Mg#,CaO/Al2O3, NaK#]
% FirstMeltPb = pressure of the first melt
% FirstAt2percent = pressure of the first extracted melt
% FirstMeltIndex = index of the first melt
% FirstAt2Index = index of the first extracted melt
% densities = [averageliquiddensity(1) averagesoliddensity averageresiduedensity averageliquiddensity(2) average_pooledliquid_density_novolumechangePT];
% LastMeltPb = Pressure of the last melt
% LastMeltIndex = Index of the last melt
% crustalthickness = crustalthickness calcualted for a column melt using density
%                   DEPTH.*initialdensity./average_pooledliquid_density_novolumechangePT./(1+(1-max(ExtFrac))./max(ExtFrac));

% poolCL_Iv = pooled trace element compostiions of melts weighted by melt velocity
% poolCL = pooled trace element compostiions of melts not weighted by melt velocity
% poolBAI4_Iv  = pooled extracted melts weighted by melt velocity
% poolBAI4  = pooled extracteds notweighted by melt velocity
% running_mBAI4 = pooling extracted melts weighted by melt velocity
% running_CL_OUT = pooling trace element compostiions of extracted melts weighted by melt velocity


if any(FirstAt2percent) ==0
    FirstAt2percent = 0;
end
if any(FirstMeltIndex) ==0
    FirstMeltIndex = 0;
end

if any(FirstMeltPb) ==0
    FirstMeltPb = 0;
end

if any(MeltExtractedPb) ==0
    MeltExtractedPb = 0;
end

if any(LastMeltIndex) ==0
    LastMeltIndex = 0;
end

if any(LastMeltPb) ==0
    LastMeltPb = 0;
end
if any(aluminumOutFrac) ==0
    aluminumOutFrac = 0;
end

if any(cpxOutFrac) ==0
    cpxOutFrac = 0;
end
if any(FirstAt2Index) ==0
    FirstAt2Index = 0;
end

if any(crustalthickness) ==0
    crustalthickness = 0;
end

if any(aluminumOutIndex) ==0
    aluminumOutIndex = NaN;
    aluminumOutPb =    NaN;
    aluminumOutFrac = NaN;
else
    aluminumOutPb =    Pb(aluminumOutIndex);
end


if any(cpxOutIndex) ==0
    cpxOutIndex = NaN;
    cpxOutPb =    NaN;
    cpxOutFrac = NaN;
else
    cpxOutPb =    Pb(cpxOutIndex);
end

phasesOutInfo = [aluminumOutPb aluminumOutFrac cpxOutPb cpxOutFrac];

assignin('caller', sprintf('ExtFrac_%g', runnumber), ExtFrac)
assignin('caller', sprintf('aluminumOutIndex_%g', runnumber), aluminumOutIndex)
assignin('caller', sprintf('cpxOutIndex_%g', runnumber), cpxOutIndex)
assignin('caller', sprintf('TransitionSp2Plag_index_%g', runnumber), TransitionSp2Plag)
assignin('caller', sprintf('TransitionGar2Sp_index_%g', runnumber), TransitionGar2Sp)
assignin('caller', sprintf('aluminumOutFrac_%g', runnumber), aluminumOutFrac)
assignin('caller', sprintf('cpxOutFrac_%g', runnumber), cpxOutFrac)
assignin('caller', sprintf('MeltExtractedPb_%g', runnumber), MeltExtractedPb)
assignin('caller', sprintf('BAI3_%g', runnumber), BAI3)
assignin('caller', sprintf('XOW_OUT_%g', runnumber), XOW_OUT)
assignin('caller', sprintf('rXOW_%g', runnumber), rXOW)
assignin('caller', sprintf('initialBulkComp_%g', runnumber), initialBulkComp)
assignin('caller', sprintf('FirstMeltPb_%g', runnumber), FirstMeltPb)
assignin('caller', sprintf('FirstAt2percent_%g', runnumber), FirstAt2percent)
assignin('caller', sprintf('FirstMeltIndex_%g', runnumber), FirstMeltIndex)
assignin('caller', sprintf('FirstAt2Index_%g', runnumber), FirstAt2Index)
assignin('caller', sprintf('LastMeltPb_%g', runnumber), LastMeltPb)
assignin('caller', sprintf('LastMeltIndex_%g', runnumber), LastMeltIndex)
assignin('caller', sprintf('CrustalThickness_%g', runnumber), crustalthickness)
assignin('caller', sprintf('Running_mBAI4_%g', runnumber), running_mBAI4)
assignin('caller', sprintf('Running_CL_OUT_%g', runnumber), running_CL_OUT)
assignin('caller', sprintf('poolCL_Iv_%g', runnumber), poolCL_Iv)
assignin('caller', sprintf('poolCL_%g', runnumber), poolCL)
assignin('caller', sprintf('poolBAI4_Iv_%g', runnumber), poolBAI4_Iv)
assignin('caller', sprintf('poolBAI4_%g', runnumber), poolBAI4)


pool_ME =  [poolBAI4_Iv poolCL_Iv];


dx=1;
dz =  154.6073 ;
M = Lnexp;
% Calculate melt productivity and crustal thickness
Mr = Iv.*(M'/dz);  % Melt production rate c(1/s)
Mr(Mr<0) = 0;      % Eliminate values where melt productivity is negative
tMx = sum(Mr).*dz;  % total melt at each X-value
%ct = sum(tMx).*dx./(10*(2.*U(n)/secyr)); % crustal thickness changed 7/7/17 to divide by full spreading rate, not half spreading rate
%ct_ininc(k+1) = sum(sum(Mr(:,k+1:end-k).*dz)).*dx./(2.*10*U(n)/secyr);

%saves the stability field of melting
for kkzz=1:Gar2SpTransition_index
    MeltType(kkzz)=1;
end
for kkzz=Gar2SpTransition_index+1:Sp2PlagTransition_index
    MeltType(kkzz)=2;
end
for kkzz=Sp2PlagTransition_index+1:size(Pb,2)
    MeltType(kkzz)=3;
end



garnetpossible = 1:Gar2SpTransition_index;
garnetexist= find(XM_OUT(garnetpossible,6)>0);
Garnet_ct_real= sum(sum(Mr(garnetexist).*dz)); %.*dx./(2.*10*U(n)/secyr);

Garnet_ct= sum(sum(Mr(1:Gar2SpTransition_index).*dz)); %.*dx./(2.*10*U(n)/secyr);
Spinel_ct = sum(sum(Mr(Gar2SpTransition_index+1:Sp2PlagTransition_index).*dz)); %.*dx./(2.*10*U(n)/secyr);
Plag_ct = sum(sum(Mr(Sp2PlagTransition_index+1:end).*dz)); %.*dx./(2.*10*U(n)/secyr);

percentGarSpPlag_temp= [Garnet_ct ; Spinel_ct; Plag_ct];
percentGarSpPlag_temp = percentGarSpPlag_temp./sum(percentGarSpPlag_temp);
percentGarSpPlag = percentGarSpPlag_temp;

average_LnExp =  sum(sum(Mr.*ExtFrac'))/sum(Mr(:));

average_Pb= nansum(nansum(Mr.*Pb'))/nansum(Mr(:));


Garnet_proportion = percentGarSpPlag(1);
Spinel_proportion = percentGarSpPlag(2);
Plag_proportion = percentGarSpPlag(3);

if XM_OUT(1,6) ==0
    'no initial garnet'
    Garnet_proportion = 0;
    Spinel_proportion = 0;
    Plag_proportion = 0;
else if XM_OUT(Gar2SpTransition_index-1,6) ==0
        'garnet is exhausted'
        Garnet_proportion = Garnet_ct_real./sum([Garnet_ct ; Spinel_ct; Plag_ct]);
        Spinel_proportion = 0;
        Plag_proportion = 0;
    end
end
for i = 1:9
    pool_garnetME(i) = sum(sum(Mr(garnetexist,:).*BAI4(garnetexist,i)))/sum(sum(Mr(garnetexist,:)));
end

for i = 1:size(PetrogenTraceElement_Strings,2)
    pool_garnetCL(i) = nansum(nansum(Mr(garnetexist,:).*CL_OUT(garnetexist,i)))/nansum(nansum(Mr(garnetexist,:)));
end

pool_garnet = [pool_garnetME pool_garnetCL];

if FirstAt2Index == 0
    firstmelt = [BAI4(1,:) CL_OUT(1,:)];
else
    firstmelt = [BAI4(FirstAt2Index,:) CL_OUT(FirstAt2Index,:)];
end







% outputMeltMatrix uses on-axis things
outputMeltMatrix(runnumber,:)=[dFdPwaterchange porosity initialwatercontent TP afewfunrnuns_U_vary(afewfunrnuns_U)  afewfunrnuns_BulkCompositions...
    FirstMeltPb MeltExtractedPb   LastMeltPb Gar2SpTransition Sp2PlagTransition Garnet_proportion Spinel_proportion Plag_proportion crustalthickness  sum(ln) average_LnExp Tdiff ...
    phasesOutInfo...
    AdiabaticMeltingInfo...
    pool_ME];

outputGarnetMeltMatrix(runnumber,:)=[dFdPwaterchange porosity initialwatercontent TP afewfunrnuns_U_vary(afewfunrnuns_U)  afewfunrnuns_BulkCompositions...
    FirstMeltPb MeltExtractedPb   LastMeltPb Gar2SpTransition Sp2PlagTransition Garnet_proportion Spinel_proportion Plag_proportion crustalthickness sum(ln) max(ExtFrac) Tdiff...
    phasesOutInfo...
    AdiabaticMeltingInfo...
    pool_garnet];

outputMeltMatrix_FullPool(runnumber,:)=[dFdPwaterchange porosity initialwatercontent TP afewfunrnuns_U_vary(afewfunrnuns_U)  afewfunrnuns_BulkCompositions...
    FirstMeltPb MeltExtractedPb   LastMeltPb Gar2SpTransition Sp2PlagTransition percentCrust(1,1) percentCrust(2,1) percentCrust(3,1) ct  sum(ln) FP_Full Tdiff ...
    phasesOutInfo...
    AdiabaticMeltingInfo...
    pool_full poolCL_full];

outputFirstMelt(runnumber,:)=[dFdPwaterchange porosity initialwatercontent TP afewfunrnuns_U_vary(afewfunrnuns_U)  afewfunrnuns_BulkCompositions...
    FirstMeltPb MeltExtractedPb   LastMeltPb Gar2SpTransition Sp2PlagTransition percentCrust(1,1) percentCrust(2,1) percentCrust(3,1) ct  sum(ln) max(ExtFrac) Tdiff ...
    phasesOutInfo...
    AdiabaticMeltingInfo...
    firstmelt];





%         ZeroIngrowth_Pool_Th230_ininc
%         ZeroIngrowth_Pool_Ra226_ininc
%         CompleteIngrowth_Pool_Th230_ininc
%         CompleteIngrowth_Pool_Ra226_ininc

FullPool_DisEQ = [ZeroIngrowth_Pool_Th230_full ZeroIngrowth_Pool_Ra226_full CompleteIngrowth_Pool_Th230_full CompleteIngrowth_Pool_Ra226_full];
NarrowPool_DisEQ = [ZeroIngrowth_Pool_Th230_ininc(101) ZeroIngrowth_Pool_Ra226_ininc(101) CompleteIngrowth_Pool_Th230_ininc(101) CompleteIngrowth_Pool_Ra226_ininc(101)];
%DISEQ_meltmass   ;
maxDISEQ_NF = [max(ZeroIngrowth_NF_Th230) max(ZeroIngrowth_NF_Ra226) max(CompleteIngrowth_NF_Th230) max(CompleteIngrowth_NF_Ra226)];

%Important degrees of melting
Fmax_OnAxis =sum(Lnexp); % OR 1-(Wnr(end))
Fmax_OnAxis_NotExp =sum(ln); % OR Fc OR 1-(WnS(end))
Fmax_AcrossAxis=sum(Lnexp_IT);
Fmax_Min = min(Fmax_AcrossAxis(Fmax_AcrossAxis>0));


FP_OnAxis = variablypooledLnExp(101); % OR Lnexp_cumulative_IT_ininc(101);
FP_OnAxis_meltmass = ExtFrac*(Lnexp./sum(Lnexp))';

%      FP_Full  : fully pooled Fmean using melt productivity
%      Lnexp_cumulative_IT_ininc(101) : on-axis pooled Fmean using melt productivity
%      Fmean_1D_meltmass : on-axis pooled Fmean using melt mass

%       NOTE: Lnexp_cumulative_IT_ininc should be the same as variablypooledLnExp!
%               Both are the Fmean up to a certain distance off-axis
% 1 = Full Pool; 101 = On-Axis

FullPoolInfo =[ct2_ininc(1)  percentCrust(1,1) percentCrust(2,1) percentCrust(3,1)];

NarrowPoolInfo = [ct2_ininc(end) Garnet_proportion Spinel_proportion Plag_proportion  variablypooledLnExp(101)];
NarrowPoolInfov2 = [Garnet_proportion Spinel_proportion Plag_proportion];



F_INFO = [FP_Full FP_OnAxis Fmax_OnAxis Fmax_OnAxis_NotExp];

% info for all runs:
StandardInfo = [...
    dFdPwaterchange porosity initialwatercontent TP afewfunrnuns_U_vary(afewfunrnuns_U)  afewfunrnuns_BulkCompositions...
    FirstMeltPb MeltExtractedPb LastMeltPb Gar2SpTransition Sp2PlagTransition phasesOutInfo Tdiff crustalthickness...
    AdiabaticMeltingInfo...
    FullPoolInfo ...
    NarrowPoolInfov2...
    F_INFO...
    FullPool_DisEQ maxDISEQ_NF...
    ];

% size(Gar2SpTransition)
% size(Sp2PlagTransition)
% size(phasesOutInfo)
% size(Tdiff)
% size(crustalthickness)
% size(AdiabaticMeltingInfo)
% size(FullPoolInfo)
% size(NarrowPoolInfov2)
% size(F_INFO)
% size(FullPool_DisEQ)
% size(maxDISEQ_NF)
% size(StandardInfo)




% TableTest = table(dFdPwaterchange,porosity,initialwatercontent,TP ,afewfunrnuns_U_vary(afewfunrnuns_U),  afewfunrnuns_BulkCompositions,...
%     FirstMeltPb, MeltExtractedPb ,LastMeltPb, Gar2SpTransition, Sp2PlagTransition, phasesOutInfo, Tdiff, crustalthickness,...
%     AdiabaticMeltingInfo,...
%     FullPoolInfo, ...
%     F_INFO,...
%     FullPool_DisEQ, maxDISEQ_NF);


%on axis: sum(ln) average_LnExp crustalthickness


%Variablly pooled melts
PoolInfo =  [distance_ininc' Lnexp_cumulative_IT_ininc percentCrust' allThickness'];
percentpools = [1:-1*1./(size(columnsWithMelt,2)-1):0];
StandardInfoRep_Pools= repmat(StandardInfo,size(PoolInfo(columnsWithMelt,:),1),1);

MASTER_POOLS_INFO =[MASTER_POOLS_INFO;[StandardInfoRep_Pools percentpools' PoolInfo(columnsWithMelt,:)]];
MASTER_POOLS_MAJORS =[MASTER_POOLS_MAJORS; poolME_ininc(columnsWithMelt,:)];
MASTER_POOLS_TRACE =[MASTER_POOLS_TRACE; poolTE_ininc(columnsWithMelt,:)];
temp_Pool = [ZeroIngrowth_Pool_Th230_ininc; ZeroIngrowth_Pool_Ra226_ininc; CompleteIngrowth_Pool_Th230_ininc; CompleteIngrowth_Pool_Ra226_ininc]';
MASTER_POOLS_DISEQ =[MASTER_POOLS_DISEQ; temp_Pool(columnsWithMelt,:)];


%Column pooled melts
%                          Pool_column_TE
%                         Pool_column_ME
%                          variablypooledLnExp
%                     variablypooledPb
%                     variablypooledTs
ColumnNum = 1:size(variablypooledLnExp,2);
ColumnPoolInfo =  [ColumnNum' variablypooledLnExp' variablypooledPb' variablypooledTs'];
%percentpools = [1:-1*1./(size(columnsWithMelt,2)-1):0];
StandardInfoRep_ColPools= repmat(StandardInfo,size(ColumnPoolInfo(columnsWithMelt_ALL,:),1),1);
MASTER_column_POOLS_INFO =[MASTER_column_POOLS_INFO;[StandardInfoRep_ColPools ColumnPoolInfo(columnsWithMelt_ALL,:)]];
MASTER_column_POOLS_MAJORS =[MASTER_column_POOLS_MAJORS; Pool_column_ME(columnsWithMelt_ALL,:)];
MASTER_column_POOLS_TRACE =[MASTER_column_POOLS_TRACE; Pool_column_TE(columnsWithMelt_ALL,:)];
%FIX: SKB
%                         temp_column_Pool = [ZeroIngrowth_Pool_Th230_ininc; ZeroIngrowth_Pool_Ra226_ininc; CompleteIngrowth_Pool_Th230_ininc; CompleteIngrowth_Pool_Ra226_ininc]';
%                         MASTER_column_POOLS_DISEQ =[MASTER_Column_POOLS_DISEQ; temp_column_Pool(columnsWithMelt,:)];
%

%Near-fractional info
if porosity == 1
    Lnexp = ln;
end
for kl = 1:size(Lnexp,2)
    Lnexp_cumulative(kl) = sum(Lnexp(1:kl)); % using (1-Wnr) has floating point issues, so use this instead
    Ln_cumulative(kl) = sum(ln(1:kl)); % using (1-Wnr) has floating point issues, so use this instead
end
Lnexp_cumulative(Lnexp==0)=0;
Ln_cumulative(Lnexp==0)=0;
indiciesWithMelt = find(Lnexp>0);

Pb_melts = Pb(indiciesWithMelt);
IncrementalMeltInfo= [Pb' T' Ts' Ln_cumulative' Lnexp_cumulative' MeltType'];
StandardInfoRep= repmat(StandardInfo,size(IncrementalMeltInfo(indiciesWithMelt,:),1),1);

MASTER_INSTANTS_MAJORS =[MASTER_INSTANTS_MAJORS; BAI4(indiciesWithMelt,:)];
MASTER_INSTANTS_TRACE =[MASTER_INSTANTS_TRACE; CL_OUT(indiciesWithMelt,:)];
MASTER_INSTANTS_INFO =[MASTER_INSTANTS_INFO; [StandardInfoRep IncrementalMeltInfo(indiciesWithMelt,:)]];
temp_NF_DISEQ=[(ZeroIngrowth_NF_Th230) (ZeroIngrowth_NF_Ra226) (CompleteIngrowth_NF_Th230) (CompleteIngrowth_NF_Ra226)];
MASTER_INSTANTS_DISEQ =[MASTER_INSTANTS_DISEQ; temp_NF_DISEQ(indiciesWithMelt,:)];

%%

MASTER_Residue_TRACE =[MASTER_Residue_TRACE; CR_OUT(indiciesWithMelt,:)];
MASTER_SolidResidue_TRACE =[MASTER_SolidResidue_TRACE; CS_OUT(indiciesWithMelt,:)];
MASTER_CPX =[MASTER_CPX; CS_minerals_OUT(indiciesWithMelt,:)];
MASTER_INSTANTS_Ds =[MASTER_INSTANTS_Ds; Ds2save(indiciesWithMelt,:)];


%%


pool_garnet_SAVE = pool_garnet;
extract_SM_1D
testXX(runnumber,:)=size(importantCompositions);
importantCompositions_Majors  = importantCompositions(1:11,1:9);
importantCompositions_Trace  = importantCompositions(1:11,10:end);
importantProportions= [...
    1 0 0; 1 0 0; 1 0 0
    0 1 0; 0 1 0; 0 1 0
    0 0 1;  0 0 1;  0 0 1
    FullPoolInfo(2:end)
    NarrowPoolInfov2];




MASTER_SM_MAJORS =[MASTER_SM_MAJORS;importantCompositions_Majors];
MASTER_SM_TRACE =[MASTER_SM_TRACE;importantCompositions_Trace];
MASTER_SM_INFO =[MASTER_SM_INFO;repmat(StandardInfo,11,1) importantFN importantProportions];



MASTER_ALL_RESULTS(runnumber,:) = StandardInfo;
MASTER_SM_DISEQ = [MASTER_SM_DISEQ; importantDISEQ];


%%


% MASTER_POOLS_INFO
% MASTER_POOLS_MAJORS
% MASTER_POOLS_TRACE
% MASTER_INSTANTS_MAJORS
% MASTER_INSTANTS_TRACE
% MASTER_INSTANTS_INFO
%

if mod(runnumber,31)==0
    eval(sprintf('save(''%s'',''outputMeltMatrix'',''outputGarnetMeltMatrix'',''outputMeltMatrix_FullPool'',''outputFirstMelt'',''MASTER_SM_MAJORS'',''MASTER_SM_TRACE'',''MASTER_SM_INFO'',''MASTER_POOLS_INFO'', ''MASTER_POOLS_DISEQ'', ''MASTER_INSTANTS_DISEQ'',''MASTER_POOLS_MAJORS'', ''MASTER_POOLS_TRACE'', ''MASTER_INSTANTS_MAJORS'', ''MASTER_INSTANTS_TRACE'', ''MASTER_INSTANTS_INFO'',''MASTER_Residue_TRACE'',''MASTER_SolidResidue_TRACE'',''MASTER_CPX'',''MASTER_INSTANTS_Ds'')',[pwd '/' subfolder '/' worksheetName2Save '.mat']))
end




assignin('caller', sprintf('AveragePb_%g', runnumber), average_Pb)      ;
assignin('caller', sprintf('AllAveragePb_%g', runnumber), variablypooledPb) ;
assignin('caller', sprintf('PooledPb_%g', runnumber), Pb_Melting_ininc) ;
assignin('caller', sprintf('PooledT_%g', runnumber), T_Melting_ininc) ;



assignin('caller', sprintf('poolBAI4_Garnet_%g', runnumber), pool_garnet)
assignin('caller', sprintf('garnetproportion_%g', runnumber), Garnet_proportion)
assignin('caller', sprintf('spinelproportion_%g', runnumber), Spinel_proportion)
assignin('caller', sprintf('plagproportion_%g', runnumber), Plag_proportion)



%mat2clip([outputMeltMatrix outputFirstMelt])

%FileNameOut = sprintf('Water_Azores_Rfrac%sWater%dTmp%d%s_dFdPchamge',rFrac_crit_str,initialwatercontent,TP,Tbchar,dFdPwaterchange)
%    Garnet_proportion
%    MeltExtractedPb
%    FirstMeltPb
%

ColumnPooled = [min(variablypooledLnExp) max(variablypooledLnExp) ...
    min(variablypooledTs) max(variablypooledTs) ...
    min(variablypooledPb) max(variablypooledPb)...
    min(variablypooledTs)-1.5.*min(variablypooledPb)...
    max(variablypooledTs)-1.5.*max(variablypooledPb)];


TpMeltEnd = [Ts(dq(bq)) - 1.5*LastAdiabaticMeltPressure];
TpMeltEnd_byF = TP-F_INFO(3).*600;


ForRevPet =[TP TpMeltEnd_byF...
    F_INFO(3)...
    LastMeltPb MeltExtractedPb   ...
    Ts(Pb==LastMeltPb) Ts(Pb==MeltExtractedPb)...
    Ts(Pb==LastMeltPb)-1.5.*LastMeltPb...
    Ts(Pb==MeltExtractedPb)-1.5.*MeltExtractedPb...
    F_INFO(1:2)...
    min(Pb_Melting_ininc) max(Pb_Melting_ininc)  ...
    min(T_Melting_ininc) max(T_Melting_ininc)  ...
    min(T_Melting_ininc)-1.5.*min(Pb_Melting_ininc)...
    max(T_Melting_ininc)-1.5.*max(Pb_Melting_ininc)...
    ColumnPooled...
    TpMeltEnd];



